\documentclass[12pt,fleqn]{article}
\title{EvoDevo Walkthrough}
\author{Willaim Fortin, Matt Helgen, Nilufer Onder, Kara Vogel}
\usepackage{times}
\usepackage{fullpage}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{algorithmic}
\usepackage{algorithm}
\usepackage{wrapfig}
\usepackage{setspace}
\usepackage{tocloft}
\usepackage{subfig}



\setlength{\parindent}{0in}
\setlength{\parskip}{0.1in}

\begin{document}

\thispagestyle{empty}
\maketitle
\newpage

\begin{onehalfspacing}
\tableofcontents
\end{onehalfspacing}

\newpage
\section{Introduction}

An important principle in biochemistry is the rise of highly complex behavior from interactions of simple rules.
However, even simple models of biochemical evolution prove challenging due to the scale of possible interactions \cite{mason:707}. 
The concentration of a particular interesting protein might be dependent on the presence and concentration of several other proteins, each of which have similar interdependencies.
By modeling the effects of these interactions, an automated process can produce configurations of a biochemical system with predefined desirable properties \cite{francoishakim}. 
This program represents a prototype to eventually produce biologically inspired situations such as oscillations or bistable switches. See section \ref{sec:future} for planned future work. 

Section \ref{sec:gettingstarted} explains how to download and compile the code.
In Section \ref{sec:framework}, we describe the basic workings of the simulation at runtime.
Section \ref{sec:interactions-mutation} details the biologically inspired mutations that may occur in the Cell.
Section \ref{sec:interactions-implementation-complex} explains the construction and function of Protein-Protein Complex interactions.
Our method and algorithm for numerical integration, Runge-Kutta, is shown in Section \ref{sec:rungekutta}. 
The evaluation function used to rank cells after rounds of numerical integration is defined in Section \ref{sec:scoring}.
The output produced by \tt EvoDevo \rm is discussed in Section \ref{sec:output}.
Current bugs in Version 1.0 and planned future work for Version 2.0 are listed in Section \ref{sec:future}.

In silico evolution programs in the past have contributed to  biological understanding through algorithmic computation with various points of focus: mutations, resource competition, robustness, protein bistable switches, and oscillations \cite{Lenski1999}\cite{Adami2006}\cite{Neyfakh2006}\cite{Bergman2003}\cite{Wagner2007} to name a few.
Integrating and expanding simple synthetic models to evolve new biological machinery such as the Goodwin model or Gillepsie stochastic interactions or the circadian clock in humans, on a computational level, contributes to sectionalized understanding of larger networks.  
Reliable small modular algorithms may one day advance and integrate to a point where certain underlying genomic attributes can be recognized at various points of potential for evolutionary adaptation or for precision gene therapy. 
 
Consider as one specific example of potential the powerful example and application of bistable switch architecture which allowed an engineered bacterial strain to invade cancer cells under hypoxic conditions in a tumor \cite{Agapakis2009}.  
Mathematical models, similar to biological examples of synthetic networks, improve the depth of understanding and provide powerful experimental tools.  
In a recent report by Buchler and Cross, sequestration of a basic leucine zipper transcription factor by a dominant-negative inhibitor confers a graded transcriptional response to a sharply ultrasensitive response in a synthetic gene network in budding yeast \cite{Buckler2009}.
The synthetic network, though engineered, facilitated understanding of sequestration, just as synthetic gene networks facilitate understanding of underlying biological gene and protein interactions, etc., and modular evolution allows rapid experimentation and exploration of potential evolutionary pathways and parameters such as the rates involved in generating stable interactions-here oscillations.  Endless possible networks, and approaches to modeling evolution have been and remain to be explored.  


\section{Getting Started}
\label{sec:gettingstarted}
The \tt EvoDevo \rm program is written in C++ to be compiled on Linux / MacOS. 
It can also be compiled on Windows using the Cygwin POSIX emulation layer.
To download the most recent source code and documentation using svn, run these commands from the command line:
\tt
\newline \newline
svn checkout http://evolvedevo.googlecode.com/svn/trunk/ evolvedevo
\newline
cd evolvedevo/EvoDevo/src
\newline
make
\newline
./EvoDevo
\rm
\newline \newline
The output folder is created in the current directory. 
Each run of the code produces a new subdirectory in the output directory named by the current process id.
Each subdirectory contains a score log, and subdirectories for each cell. 
The score.log file details the best scored cell at each generation. 
Within the cell subdirectories, an output file is produced for each generation. 
This file shows the configuration of the cell at that generation, including all of the kinetic rates and interactions. 
It also contains the data produced by Runge-Kutta. 
For generations during which the Cell is scored as the best cell, it also has a seperate output file with the precise Runge-Kutta output, and several graphs based on that output.



\section{Simulation Framework}
\label{sec:framework}


The kinetics of each Cell can be represented by a graph structure. 
The initiaization for this example is presented in Fig. \ref{fig:graphexample}. 
Molecules are represented by nodes, and edges represent interaction and relationships between molecules.  

The rate of transcription for $A$ ($k_{txn\_A}$) is represented by the edge between $gA$ and $mA$. 
Likewise, the rate of translation $k_{tsln\_A}$ is the edge between $mA$ and $A$.  
In our cellular model, we assume the nucleic acid and molecular machinery necessary for transcription and translation is available in abundance. 

In addition, the $\emptyset$ node is not a molecule, but represents the spontaneous degradation and utilization of molecules, which takes place at a rate represented by $k_{deg\_(name)}$.

Consider a simulation consisting of 100 sample cells run for1000 generations, with scoring every ten generations. For the program, this means several things:

\begin{figure}[b!]
\centering
\includegraphics[scale=0.35]{figures/CellDiagram2.png}
\caption{Graph representation of an example cell}
\label{fig:graphexample}
\end{figure}


\begin{enumerate}
\item100 independent \tt Cell\rm  objects are created (\it numCells \rm = 100)
\item Each Cell will can be initialized with two \tt Protein \rm objects $A$ and $B$, two \tt mRNA \rm objects $mA$ and $mB$, and two \tt DNA \rm objects $gA$ and $gB$ (\it minBasic \rm = 2)
\item Basal rates of transcription and translation are established for each group ($gA, mA, A$)
\item Each Cell will undergo 1000 rounds of possible mutation, and will be scored after every ten rounds (\it scoringInterval = 10\rm, \it maxGenerations=1000 \rm = 6)
\item Scoring will be assess properties of the ODE's and the nearness to biological switch, oscillatory or bifurcation behavior.
\item Data accumulated after each (\it scoringInterval \rm = 2)l will be output to a folder in a series of files for each bestCell.
\end{enumerate}



\section{Class Breakdown}
\label{sec:classstructure}
See the class diagram toward the end of this document to supplement this section.

The major overarching class in EvoDevo is the \tt Experiment \rm class. It is responsible for storing all the cells, and all of the information about the experiment as a whole, as well as for conducting the experiment. It stores all the \tt Cell \rm objects in a vector, and a few other data (such as the number of elapsed generations) are stored in member variables.

An \tt Experiment \rm holds a number of \tt Cell\rm  objects, the basic unit of the experiment. The \tt Cell \rm holds all the \tt Interaction \rm and \tt Molecule \rm objects for the experiment, and manages all its internal operations (calculations derived from its properties, mutation of new properties, and management of its data).

The \tt Molecule \rm objects are responsible for handling all long-term and temporary values associated with their respective molecule. This includes names and ID numbers used both for distinction by the program and for reference by the user, and concentration values at various stages of algorithmic computation.

The actual \tt Molecule \rm objects used in the program are actually subclasses of Molecule. These subclasses are \tt Protein\rm , \tt MRNA\rm , and \tt DNA\rm . Of the three, Protein is the only one that currently requires special properties (these being Post-Translational Modifications and a subunit count), but some possible future functions will require different classes for each.

The \tt Interaction \rm class's purpose is to define the rate at which a molecule produces another molecule, decomposes into a different molecule, or spontaneously degrades. An \tt Interaction \rm holds three important variables: rate, productFactor, and substrateFactor. 
The productFactor and substrateFactor are a constant multiplication of the amount of concentration added to the product \tt Molecule \rm or removed from the substrate \tt Molecule \rm , and are initialized to 1.0. 
For example, mRNA is not consumed when producing Proteins, so the substrateFactor for that interaction would be set to 0.

The Molecule and Interaction classes are related by a directed graph structure, with \tt Molecule \rm nodes, and \tt Interaction \rm edges. The graph provides structure, and directly relates each \tt Interaction \rm to the \tt Molecule \rm objects which are interacting. Most interactions are dependent only on a single substrate, and a single product, but some interactions might need to take two substrates into account. See Section \ref{sec:interactions-implementation-complex} for more on Protein-Protein Complex interactions.

The \tt Quantity \rm object holds a pointer to a float, and a character array allowing a float value to be given a name or tag.

A \tt Term \rm object holds a vector of \tt Quantity \rm objects, and provides a method \tt Term::getValue() \rm which returns the product of it's quantity. This is a useful construct for setting up complex interactions as described in Section \ref{sec:interactions-implementation-complex}.



\section{Interactions}
\label{sec:interactions}
Simple systems initialize with minimal user defined basic proteins and "evolve" by mutations which increase interaction and cell components (repressors, complexation)

\subsection{Evolution of Cells: Mutation Types}
\label{sec:interactions-mutation}

\begin{figure}[h]
\centering
\includegraphics[scale=0.35]{figures/MolecularInteractions.jpg}
\caption{Molecular Interactions Schematic}
\label{fig:interactions}
\end{figure}

During each generation, each cell will undergo a mutation chosen by \tt Cell.mutate()\rm. 
This mutation will add or change values to the ODE's describing the cell. A random number \tt category \rm is drawn first to choose a category for possible mutations, and a second random number \tt subcategory \rm is drawn to choose the specific mutation type.

During each generation, each cell will undergo a mutation chosen by \tt Cell.mutate()\rm. This mutation will add or change values to the ODE's describing the cell.


There are three categories of possible mutations and nine mutation types:
\renewcommand{\theenumi}{\Roman{enumi}}
\renewcommand{\labelenumi}{\theenumi.}

\renewcommand{\theenumii}{\arabic{enumii}}
\renewcommand{\labelenumii}{\theenumii.}

\begin{enumerate}
	\item Small Mutation
	\begin{enumerate}
		\item Forward Rate Change - \tt Cell.fwdRateChange() \rm - \it A randomly chosen interaction representing a production, translation, or transcription rate is assigned a new random \tt value\rm.\it A forward rate involves the production of a molecule, such as translation, transcription, or production of a PTM or Complex Protein. \rm
		\item Reverse Rate Change - \tt Cell.revRateChange() \rm - \it A randomly chosen interation representing a reverse rate is assigned a new random \tt value\rm. \it A reverse rate involves a PTM or Complex protein disassociating to it's original molecules.\rm
		\item Degredation Rate Change - \tt Cell.degRateChange() \rm - \it A randomly chosen interaction representing a degradation rate is assigned a new random \tt value\rm. \it A degradation rate involves any molecule which is being degraded and is lost.\rm
		\item Histone Modification - \tt Cell.histoneMod() \rm - \it A randomy chosen interaction representing a translation rate is assigned a new random \tt productFactor \rm \it between [0,2]. \rm
		\item Post-Translational Modification - \tt Cell.newPostTransMod() \rm - \it A new protein is created as a production of a previously existing randomly chosen protein. The new protein receives a copy of the \tt postTransMods \rm \it array, with a random index incremented by 1, representing the new modification.\rm
	\end{enumerate}
	
	\item Large Mutation
	\begin{enumerate}
	\setcounter{enumii}{6}
		\item New Gene/Basic Protein - \tt Cell.newBasic() \rm - \it A new DNA, mRNA, and Protein are added to the cell with randomly set kinetic rates. \rm
		\item Protein-Protein Complex - \tt Cell.newComplex() \rm - \it Two randomly selected proteins are selected to complex together to form a new protein. The Protein-Protein Complex which is formed depends on the concentration of each of it's component Proteins. See Section \ref{sec:interactions-implementation-complex} for a detailed walkthrough. \rm
		\item Protein-Promoter Complex - \tt Cell.newPromoter() \rm - \it Protein-Promoter Complexes are not implemented in \tt EvoDevo v1.0 \rm.
	\end{enumerate}
	\item Null Mutation - \it A null mutation does not affect the \tt Cell \rm in any way, and could occur naturally, or as a result of a failed mutation.\rm
\end{enumerate}

A Small Mutation has smaller effect to the system and can be implemented with a higher probability than a Large Mutation.

When a mutation is being selected, one of the categories will be selected randomly, with a probability based on user input. 
Once a category is chosen, a more specific type of mutation from that category will be further selected randomly. 
Once that specific type is chosen, the target (e.g., which two proteins to combine in a Protein-Protein Complex) and effect (e.g. the reaction rate) of that mutation will be randomly generated and applied to the cell. 

Some mutations may be invalid, such as adding a new gene to a cell that is at the \it maxProtein \rm limit, or a Protein-Protein complex which already exists. 
In this case, the mutation can either be ignored as a Null mutation (slow evolution) or redrawn until a valid mutation is selected (fast evolution). 
This option is defined by the user.

\subsection{Mathematical Model}
\label{sec:interactions-equations}
We use a deterministic mass-action model, in which the concentration of molecules in the cell at a particular time are dependent on the concentration at the previous timestep and kinetic rates associated with an interaction \cite{edelstein1988}. 

\input{equations/equ}

\subsection{Implementation}
\label{sec:interactions-implementation}
A \tt vector<Molecule*> \rm\it proteinList \rm in the \tt Cell \rm class holds references to all proteins in the graph, for easy access to a randomly selected protein.

\subsubsection{Protein-Protein Complex}
\label{sec:interactions-implementation-complex}
The complexation of $P1$ and $P2$ into $P1P2$, at the rate of

$ \frac{d[P_1P_2]}{dt} = [P_1] \times [P_2] \times k_{complex} - [P_1P_2] \times k_{r\_complex} - [P_1P_2] \times k_{deg\_complex} $

is represented as follows: $P_1$, $P_2$, and $P_1P_2$ are \tt Protein \rm objects, and each edge is represented as an \tt Interaction \rm object.

Each \tt Protein \rm object contains a \tt Quantity \rm object \it concentration\rm, and an \tt integer \rm array \it postTransMods \rm denoting the number of Post-Translational Modifications (PTMs) it has. 
$P_1$ and $P_2$ may or may not have PTMs.

Each \tt Interaction\rm contains a \tt Term \rm object \it rate\rm, which refers to the rate at which it occurs, and two \tt float \rm values \it substrateFactor \rm and \it productFactor \rm indicating how the rate of the reaction is scaled in calculating its effect on the concentration of its substrate and its product, respectively. 
A value of 0.5 for \it productFactor \rm signifies that half of the calculated value of the rate is applied to the product.

The edges going from $P_1$ or $P_2$ to $P_1P_2$, the synthesis edges, store references to the same \tt Quantity \rm object, as their rates are governed by the same constant $k_{complex}$. The degradation edges (those going from $P_1P_2$ to $P_1$ or $P_2$ also store \tt Quantity \rm objects with identical constants. It is important to node that these \tt Quantity \rm objects are instances of the same struct used by the \tt Quantity \rm objects in the \tt Protein \rm objects to store their concentrations, even though one stores constant values and the other does not.


To create a complex, the program finds two \tt Protein \rm instances within a cell that have no common complex. 
The program assumes a preexisting complex between $X$ and $Y$ if they both have outgoing edges to another \tt Protein \rm $Z$, although it may be $X$ and $Z$ complex into $Y$, or $Y$ and $Z$ complex into $X$. 
Once \tt Protein \rm $A$ and $B$ are found, a third \tt Protein \rm $AB$ is created, and the four necessary edges, two representing degradation of $AB$ into $A$ and $B$ respectively and two representing the binding of $A$ and $B$ into $AB$ are created. Rate constants are determined for the degradation and binding reactions and stored in seperate \tt Quantity \rm objects, which are stored in \tt Term \rm objects. 
To the \tt Terms \rm for the binding rates of $A$ and $B$ a reference to the other molecule's \it concentration \rm must be added. 
That is, the \it rate \rm variable in $A$ represents $k_{complex_{AB}} \times [B]$, while \it rate \rm in $B$ represents $k_{complex_{AB}} \times [A]$. 
This is done because the Runge-Kutta algorithm already factors in each molecule's own concentration, so in this way they will both be equivalent to $k_{complex_{AB}} \times [A] \times [B]$. 
The \it rate \rm for the degradation of $AB$ only varies with $[AB]$; no data other than the constant $k_{deg_{AB}}$ must be stored in it. 
Then, in order to keep the reactions synchronized, a \it substrateFactor \rm and \it productFactor \rm must be set.

Because one $A$ and one $B$ combine to form only one $AB$, only half as much of the $A$ and $B$ that is lost can be added as $AB$, so \it productFactor \rm must be 0.5 for the binding edges. 
It would also work to set one as 1.0 and the other as 0.0, or as any other pair of numbers that summed to exactly 1.0, because the rate of reaction is equal for both. 
In the degradation edges, one $AB$ degrades to form one $A$ and one $AB$, so \it substrateFactor \rm must be 0.5, in order to avoid removing too much $AB$. 
Although it would work to keep only one of the two fctors and simply scale the rate, e.g. a \it substrateFactor \rm of 0.5 has the same effect as a halved rate and a \it productFactor \rm of 2.0, we preferred to keep the rate, which has a real-world equivalent, free of programming manipulations.


\section{Runge-Kutta Evaluation}
\label{sec:rungekutta}
To approximate the concentrations of the molecules in the cell over time, the 4th order Runge-Kutta method will be used. 
Runge-Kutta is an efficient, accurate method to approximate solutions to a set of Ordinary Differential Equations (ODEs)\cite{EdwardsPenney}. 
Although the result is identical to the standard Runge-Kutta method, we describe our Runge-Kutta algorithm, which takes advantage of the graphical structure of our biochemical system, from which node and edges formulate derivatives.

\subsection{Forming Derivatives: Using Node-Edge relationship to build a system of ODEs}
\label{sec:rungekutta-buildode}
During each Runge-Kutta iteration, all of the edges of the graph are traversed in no particular order. 
As the terms of the equations are simultaneous and commutative, it only matters that all of them are accounted before the next iteration.
Molecules whose concentration values will be considered as constant (such as genes and the "$\emptyset$" molecule) do not need to be approximated using the Runge-Kutta algorithm and will have no corresponding $rk$ values. 

For each molecule, incoming edges are considered as additive interactions, and outgoing edges as negative. i
The effect of the interaction is calculated by multiplying the concentration of the substrate molecule by the rate of he interaction.  
Each interaction represents a term in the differential equation of that molecule, and summing all of the effects results in the differential equation for that molecule.


\subsection{Solving the ODE using Runge-Kutta}
\label{sec:rungekutta-solve}

At each time $t$, Runge-Kutta calculates four intermediate values for each molecule, $(rk_1 .. rk_4)$ for use in approximating the concentration at time $t+h$. 
The first value, $rk_1$, is dependent only on the concentrations of the molecules at time $t$, and each subsequent $rk$ value is dependent on the concentration, the previous $rk$ value, and the timestep $h$ ($rk_2$ depends on $rk_1$, $rk_3$ on $rk_2$ etc).

The algorithm loops through the molecules four times, and calculates one of the intermediate $rk$ values on each pass. This is done as follows:

For each molecule, incoming edges are considered as additive interactions, and outgoing edges as negative. 
The algorithm loops through all of these edges, and adds or subtracts from the current molecule's $rk$ value. 

Consider calculating $rk_1$ for $mA$: 

$\frac{d[mA]}{dt} = 1.0 \times k_{txn\_A} \times [gA] - 1.0 \times k_{deg\_mA} \times [mA] - 0.0 \times k_{tsln\_A} \times [mA]$
\newline

Note the factor 0.0 in the final term, as mRNA is not consumed in the translation process.

To calculate $mA_{rk_1}$ from this graph, each edge corresponds to a term in $mA$'s differential equation, and must be applied: 

$mA_{rk_1} = mA_{rk_1} + 1.0 \times k_{txn\_A} \times [gA]$
\newline
$mA_{rk_1} = mA_{rk_1} - 1.0 \times k_{deg\_mA} \times [mA]$
\newline
$ma_{rk_1} = mA_{rk_1} - 0.0 \times k_{tsln\_A} \times [mA]$

The sum of these effects is equivalent to the differential equation described above.

After one pass through all of the molecules, every molecule has an $rk_1$ value, which is a basic approximation of what the change in concentration will be during the next timestep.

After a second pass through all of the molecules, every molecule has an $rk_2$ value, which is an alternate approximation of the change in concentration which takes the $rk_1$ value into account. 
This continues similarly for $rk_3$ and $rk_4$ according to the standard Runge-Kutta numerical approximation until all molecules have an $rk_4$ value.

Finally, in order to make an accurate approximation of the actual change in concentration between time $t$ and time $t+h$, we loop through all of the molecules and take a weighted average of thier $rk$ values to calculate change in concentration.

$Concentration(t+h) = Concentration(t) + h/6 \times (rk_1 + 2(rk_2) + 2(rk_3) + rk_4)$


\section{Scoring}
\label{sec:scoring}
Every $n$ rounds of mutation (\it generation \% n = 0\rm), the Cells will go through a scoring phase, where $n$ is chosen by the user.

Scoring acts as an initial screening providing a set number of cells from a large population automatically. 
The function selects cells which display nearness to desired mathematical behavior, for further analysis.

The Runge-Kutta evaluation will produce a set of points, $(t,concentration)$ for each molecule in the cell. 
These points are used as input to one of two screening functions which score the existing system for each cell based on nearness to oscillatory or bifurcation behavior (bistable switch).  
After sorting the \tt Cell \rm, these will be ranked and reevaluated by Runge-Kutta at higher precision, lower step size, and an extended period.  
The data from the desired preselected molecule type (basic proteins, PTM proteins, complexed proteins or mRNA) from the bestCell will be saved to a file, as well as graphed with Gnuplot\tt (www.gnuplot.info)\rm.

\subsection{Selection for Oscillitory Behavior}
\label{sec:scoring-oscillation}
The initially implemented scoring function \tt Cell.getScore()\rm detects oscillatory behavior in a simple way. 
After each \tt Molecule\rm in the \tt Cell\rm has been evaluated by Runge-Kutta, the scoring function counts the number of minima and maxima present in the data. 
A Cell's score is equal to \tt MAX(MoleculeScores).\rm 
This quickly allows attention to be drawn to the subset of Cells which contain an oscillating Molecule.

Assuming the simulation runs for $m$ generations, and is scored every $n$ rounds, the saved file will contain the evaluation data from the best cell at generations $0, n, 2n, ..., (m-n), m$ as well as a complete description of the cell.  
The datafile output by the program contains a description of the best cell at each scoring generation. 
This description includes the proteins in the cell, the rates of all kinetic interactions, and any user defined variables or selections. 
The concentration data for the molecule(s) output by Runge-Kutta is also included, and can optionally be plotted as a graph image. 
The data is plotting using gnuplot, and can be directed to a number of file formats.

%\subsection{Selection for bistable switches}

\section{Output}
\label{sec:output}
\subsection{Directory Structure}
\label{sec:output-directory}
Each cell has its own subdirectory to hold output files at each generation. 
A cell can have output for two different evaluation modes: standard and precise. 
Every cell will have a standard output file for all scored generations, and a precise output file for generations in which it is the highest scored cell. 

\begin{itemize}
   \item output/
   \begin{itemize}
     \item cell0/ 
     \begin{itemize}
       \item gen0.std.txt
       \item gen1.std.txt
       \item gen1.precise.txt
       \item gen1.precise.complete.png
       \item gen1.precise.mol1.png
       \item gen1.precise.mol2.png
       \item mutations.log.txt
     \end{itemize}
     \item cell1/
     \begin{itemize}
       \item gen0.std.txt
       \item gen0.precise.txt
       \item gen0.precise.complete.png
       \item gen0.precise.mol3.png
       \item gen0.precise.mol4.png
       \item mutations.log.txt
     \end{itemize}
     \item score.log.txt
   \end{itemize}
\end{itemize}

\subsection{Data Plotting}
\label{sec:output-plots}
The \tt Cell\rm with the highest score at each generation is reevaluated by Runge-Kutta at higher precision, producing more points at a smaller step size. 
These points are used by Gnuplot to produce plots to visualize the data. 
A single plot showing all molecules will be produced, as well as an individual plot for each molecule. 
For a \tt Cell\rm with 6 molecules, one graph would be produced with six lines plotted, and six graphs would be produced with a single line plotted.


 
\begin{figure}[h]
  \centering
  \subfloat[Complete]{\label{fig:completeGraph}\includegraphics[width=.4\textwidth]{figures/example_complete.png}}
  \subfloat[Individual]{\label{fig:moleculeGraph}\includegraphics[width=.4\textwidth]{figures/example_mol.png}}
  \caption{Generated Plots by \tt Gnuplot \rm}
  \label{fig:plot}
\end{figure}

\section{Future Work}
\label{sec:future}
\subsection{Version 1.0 Simplifications}
\label{sec:future-v1simp}
Protein-Promoter interactions do not function.

A complex protein must be made up of two distinct proteins.

\subsection{Version 2.0 Planned Work}
\label{sec:future-v2plans}
Implement protein-promoter interactions.

Allow complex proteins to be made up of proteins complexing with themselves.

More complex evaluation function(s).

Allow parameters to be selected from command line or a config file.




\begin{figure}[h]
\centering
\includegraphics[scale=0.9]{figures/classDiagram.pdf}
\caption{Class Diagram}
\label{fig:classDiag}
\end{figure}



\newpage
\bibliographystyle{unsrt}
\bibliography{WalkthroughRefs}

\end{document}
